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1.  Introduction 


Under  the  support  of  the  US  Air  Force  Summer  Faculty  Fellowship  Program  (SFFP), 
Xiaolin  Li,  faculty  member  of  the  Stony  Brook  University,  and  Bernard  Moore,  a  graduate 
student  of  Stony  Brook  University  carried  out  a  research  visit  to  Edwards  Air  Force  Base. 
During  this  visit,  the  visiting  faculty  member  collaborated  with  Dr.  Daniel  Reasor  and  Mr. 
Keerti  K.  Bhamidipati  in  a  joint  research  effort  on  three  problems  of  potential  Air  Force 
interests. 

Our  approach  is  through  the  Lagrangian  front  tracking  method  which  is  an  adaptive 
computational  method  to  follow  the  dynamic  evolution  of  a  discontinuous  interface.  It 
represents  the  interface  as  a  set  of  topologically  connected  marker  points  and  uses  finite 
difference  method  to  each  smooth  subdomain.  The  discontinuity,  or  the  front,  is  treated  as 
an  interior  boundary  for  the  PDF  solver  in  each  subdomain.  It  is  known  that  this  method 
can  give  high  resolution  at  subgrid  level  and  keep  sharp  corners  with  acute  angles  [2,  3]. 
Over  the  last  twenty  years,  Xiaolin  Li,  in  collaboration  with  his  Stony  Brook  colleagues,  has 
applied  this  method  to  many  scientific  and  engineering  problems  including  fluid  interface 
instabilities,  fluid  structure  interactions  and  phase  transition  problems.  It  is  our  intention 
the  explore  opportunities  of  applying  this  computational  tool  to  problems  of  interests  to  the 
Department  of  Defense  including  both  Army  and  Air  Force. 

During  this  research  visit,  we  have  performed  benchmark  tests,  collaborated  in  design 
of  numerical  algorithms  and  proposed  for  future  research  problems  in  three  areas.  The  first 
problem  is  the  study  of  supersonic  and  transonic  flow  around  the  aircraft  wings.  During 
the  visit,  we  implemented  functions  to  initialize  two  dimensional  wings  from  the  database 
website  of  UIUC.  The  visiting  graduate  student,  Bernard  Moore,  carried  out  verification 
study  of  front  tracking  code  on  supersonic  flow  around  a  cylindrical  obstacle  and  compare 
the  distance  from  the  bow  shock  to  the  surface  of  the  obstacle  with  the  analytical  solution. 
The  verification  study  extended  to  different  wings,  including  the  F-16  wing  in  supersonic  flow. 
In  subsequent  study,  we  found  discrepancy  between  the  front  tracking  code  and  solution  from 
AERO  Suite  on  transonic  flow.  We  attributed  this  to  the  inaccurate  boundary  treatment  of 
the  front  tracking  code  and  will  improve  the  front  tracking  code. 

The  second  problem  is  on  the  modeling  of  fabric  surface  using  the  spring-mass  model 
and  its  application  to  the  simulation  parachute  inflation.  One  important  accomplishment  we 
made  during  this  visit  is  the  design  and  implementation  of  algorithm  to  calculate  the  von 
Mises  stress  which  is  critical  for  parachute  deployment  and  inflation.  We  met  with  Air  Force 
experts  who  gave  us  suggestions  and  recommendations  on  the  validation  study  of  parachute. 
Under  their  suggestion,  we  obtained  experimental  data  from  Prof.  Jean  Potvin  the  drag 
force  in  skydiving  experiment  of  the  C-9  personnel  parachute.  We  carried  out  numerical 
simulation  and  compared  with  the  experimental  data.  The  comparison  showed  agreement 
on  duration  and  peak  drag  of  the  C-9  inflation,  but  also  showed  some  disagreement  in  after 
inflation  drag.  This  comparison  study  provides  vital  information  on  the  improvement  of  the 
parachute  model. 

During  the  visit,  we  also  discussed  the  potential  application  of  the  phase-transition  mod¬ 
ule  in  the  front  tracking  code  to  the  airplane  icing  test.  This  newly  developed  front  tracking 
capability  to  compute  water  freezing  and  melting,  together  with  the  fluid  advection  will 
provide  a  cost-efficient  computational  platform  for  the  study  of  aircraft  icing. 
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2.  Simulation  of  Shock  Around  Airfoil 

lovnovich  and  Raveh  [5]  studied  the  oscillatory  shock-buffet  of  the  NACA-0012  and 
NACA-64A204  airfoils  and  showed  that  these  oscillations  occur  and  diminish  over  a  range  of 
attack  angle  in  transonic  flow.  For  example,  it  was  shown  that  at  Mach  0.75,  the  shock-buffet 
instability  for  the  NACA-64A204  occurs  between  the  onset  attack  angle  of  5.1  degree  and 
offset  angle  of  8.0  degree.  It  was  observed  that  shock-buffet  onset  occurs  when  the  shock  is 
just  behind  the  location  of  maximum  curvature  where  the  surface  slope  is  close  to  zero. 

Our  Edwards  advisor  Daniel  Reasor  introduced  this  phenomenon  at  the  beginning  of  our 
visit  and  recommended  that  we  use  the  front  tracking  code  to  study  the  formation  of  shock 
buffet  which  could  help  enhance  the  understanding  of  airfoil  in  the  transonic  regime.  Under 
his  suggestion,  we  implemented  the  initialization  functions  for  aircraft  wings  in  the  front 
tracking  code  during  the  first  week  and  Bernard  Moore  carried  out  a  series  of  simulations, 
first  the  verification  of  the  code  on  supersonic  flow  and  bow  shock  formation,  and  then  the 
study  of  transonic  shock  buffet. 


Figure  1:  Verification  study  of  front  tracking  code  on  supersonic  flow  around  a  cylindrical  obstacle.  The 
plots  show  the  air  pressure  in  simulations  of  four  different  Mach  numbers. 
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Mach  Number 

Thcsoretical  Distance 

Simulated  Distance 

1.25 

7,6667 

7.667 

1.50 

3,0760 

2,985 

2.00 

1,2406 

1,370 

2.50 

1,0109 

0.996 

Figure  2:  Comparison  of  the  shortest  distance  from  the  shock  front  to  the  surface  of  the  cylinder  with 
theoretical  values. 


Figure  3:  Comparison  of  the  shortest  distance  from  the  shock  front  to  the  surface  of  the  cylinder  with 
theoretical  curve. 


On  verification,  Moore  performed  simulations  on  axially  symmetric  cylinder.  He  con¬ 
tinued  simulations  for  the  NACA-0012,  NACA-0015,  WORTMANN-FX-76MP-120,  and 
NACA-64A204  airfoils.  By  running  various  simulations  for  cylinders,  we  found  an  optimal 
grid  size  which  gave  us  results  that  closely  resemble  those  previously  presented  in  [9,  8,  1]  on 
the  shape,  curvature,  and  standoff  distance  of  the  bow  shock.  The  comparison  of  supersonic 
flow  around  cylinder  gave  excellent  agreement  between  theory  and  the  numerical  solution 
from  the  front  tracking  code.  However,  on  simulations  at  transonic  speeds,  a  converging 
solution  requires  computational  mesh  20  times  of  refined  mesh  as  in  the  supersonic  cases. 
We  attribute  this  shortcoming  of  the  front  tracking  code  due  to  the  lack  of  enforcement  of 
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Figure  4:  Pressure  plot  in  simulations  of  bow  shock  formation  for  supersonic  flow  around  different  aircraft 
wings. 


Figure  5:  Pressure  plot  in  the  simulation  of  transonic  shock  on  top  of  F-16  wing  using  the  front  tracking 
method. 

mass  conservation  at  both  the  external  bonndary  and  the  aircraft  wings.  Althongh  the  study 
did  not  produce  the  desired  result,  it  points  to  the  direction  for  improvement  of  the  front 
tracking  code  for  the  simulation  of  thin  aircraft  wings. 

Figure  1,  Figure  2,  and  Figure  3  show  the  result  of  the  verification  study  on  the  shock 
formation  in  simulations  of  supersonic  flow  around  a  cylinder.  Figure  4  shows  the  formation 
of  bow  shock  in  simulations  of  supersonic  flow  around  different  aircraft  wings.  Figure  5  is 
the  pressure  profile  in  the  simulation  of  transonic  flow  around  the  wing  of  F-16  jet  fighter. 

3.  Parachute  Modeling 

Previous  studies  on  parachute  mostly  focus  on  the  steady  drag  and  velocity  in  the  terminal 
descent  regime  and  their  dependence  on  parachute  geometry,  dimension  and  payload.  In 
contrast  to  terminal  descent,  parachute  inflation  has  a  very  short  time  duration,  typically  in 
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a  few  seconds.  However,  this  short  period  regime  of  parachute  dynamics  is  a  very  important 
phase  of  the  deceleration  system.  Any  malfunction  such  as  inversion,  barber’s  pole,  or 
jnmper-in-tow  could  have  serious  effect  on  the  fate  of  the  personnel  and  cargo  delivery. 
Parachute  inflation  is  a  complex  problem  involving  aerodynamics,  structure  dynamics  and 
elasticity.  Researchers  have  studied  parachute  inflation  with  different  methods  including 
empirical  analysis,  semi-numerical  simulation  and  experiment  [10,  11]. 


Figure  6:  The  left  plot  shows  spring  model  on  a  triangulated  mesh.  Each  vertex  point  in  the  mesh  represents 
a  mass  point  with  point  mass  m.  Each  edge  of  the  triangles  has  an  equilibrium  length  set  during  initialization 
and  the  changing  length  exerts  a  spring  force  on  the  two  neighboring  vertices  in  opposite  directions.  Gores 
are  added  as  curves  with  larger  spring  constant.  The  right  shows  the  opened  parachute  through  the  spring 
model. 


Eigure  7:  Convergence  test  of  the  string  chord  under  mesh  refinement.  Erom  left  to  right  the  total  number 
of  points  in  the  string  are  50,  100,  200,  400  respectively.  The  total  mass  of  the  string  chord  as  well  as  the 
payload  at  the  lower  end  are  kept  constant  in  the  simulations. 
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Figure  8:  Convergence  test  results  of  spring  model  on  string  chord.  In  this  test  the  string  chord  is  fixed  at 
one  end  while  the  other  end  has  a  payload  and  is  free  to  move.  The  simulations  are  on  the  sequences  with 
50,  100,  200,  400  points  respectively.  The  total  mass  M  =  Nm  is  a  constant  in  the  simulation.  The  plots 
show  the  convergence  of  total  kinetic  energy  of  the  system. 


mesh  size 

ei 

Cfc 

6p 

50  and  100 

0.00374 

0.01664 

0.03464 

100  and  200 

0.00184 

0.00821 

0.01726 

200  and  400 

0.000918 

0.00406 

0.008614 

Table  1:  Convergence  tests  of  spring  model  for  a  swing  chord.  In  the  computational  sequences,  the  total 
mass  of  the  swing  chord  is  fixed.  As  the  number  of  points  increases,  the  point  mass  is  reduced  accordingly. 
Cauchy  error  is  calculated  on  two  consecutive  mesh  sequences.  Column  e/,  e/c,  and  are  errors  of  total 
length,  total  kinetic  energy  and  total  spring  potential  energy,  respective.  The  numerical  results  show  the 
first  order  convergence  for  each  of  them. 


Quantitative  analysis  of  fabric  aerodynamic  deceleration  system  (ADS)  is  a  challenging 
problem  due  to  its  complexity.  An  accurate  numerical  and  computational  platform  of  such 
system  requires  sophisticated  technology  in  computational  fluid  dynamics,  computational 
structure  dynamics  and  fluid-structure  interaction.  We  have  approached  the  problem  through 
the  spring-mass  fabric  model.  Figure  6  is  a  schematic  plot  of  the  spring-mass  fabric  model.  In 
our  previous  papers  [7,  6],  we  introduced  the  use  of  spring-mass  system  as  a  meso-scale  model 
to  mimic  the  dynamic  motion  of  fabric  surface  and  incorporate  this  model  into  the  numerical 
modeling  and  simulation  of  parachute  system.  We  showed  that  the  system  is  a  conservative 
system  (without  external  force  and  damping)  and  the  motion  is  purely  oscillatory  in  the 
direction  tangential  to  the  fabric  surface.  We  also  proved  that  there  exists  an  upper  bound 
for  the  eigen-frequency  of  such  oscillatory  motion 


CJ  < 


where  k  is  the  spring  constant,  m  is  the  point  mass,  and  M  is  the  maximum  number  of  the 
neighbors  a  spring  vertex  point  can  have.  The  spring  model  preserves  the  property  of  a  fabric 
surface  with  stiff  elasticity  but  no  bending  energy.  It  reacts  to  compression  force  in  tangential 
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mesh  size 

ga 

Cp 

15  and  30 

0.02528 

1.91810 

1.97967 

30  and  60 

0.01507 

1.21784 

1.21772 

60  and  120 

0.00604 

0.53550 

0.52837 

Table  2:  Convergence  tests  of  spring  model  for  a  fabric  drum.  In  the  computational  sequences,  the  total 
mass  of  the  membrane  is  fixed.  As  the  number  of  points  increases,  the  point  mass  is  reduced  accordingly. 
Cauchy  error  is  calculated  on  two  consecutive  mesh  sequences.  Column  e^,  Ck,  and  Cp  are  errors  of  total 
area,  total  kinetic  energy  and  total  spring  potential  energy,  respective.  The  numerical  results  show  the  first 
order  convergence  for  each  of  them. 


Figure  9:  We  use  the  von  Mises  formula  Eq.  (3)  to  calculate  the  fabric  stress  in  the  spring  model.  The  plots 
show  the  von  Mises  stress  of  rectangular  and  triangular  membranes  when  pulled  from  corners. 

direction  with  realistic  wrinkling.  However,  questions  remain  open  as  to  whether  such  system 
is  convergent  when  the  computational  mesh  is  refined,  and  if  so,  to  what  continuum  model 
it  converges.  During  this  research  visit,  we  carried  out  the  convergence  study  and  provided 
answer  to  the  first  part  of  the  question,  that  is,  the  system  is  indeed  convergent. 

We  have  performed  a  series  of  numerical  simulations  for  the  fixed  boundary  string  in  two 
dimensions  and  membrane  in  three  dimensions  under  computational  mesh  refinement.  In 
both  2D  and  3D  simulations,  we  keep  the  spring  constant  of  the  string  (2D)  or  membrane 
(3D)  as  a  constant.  We  also  keep  the  total  mass,  which  is  the  summation  of  all  point  mass 
M  —  Nm,  as  a  constant,  here  N  is  the  total  number  of  points  and  m  is  the  mass  of  each 
mass  point. 

The  Cauchy  sequences  are  numerical  solutions  with  increasing  number  of  mass  points 
(and  segment).  The  2D  sequences  are  for  total  number  of  points  N  =  50,100,200,400 
respectively.  We  recorded  characteristic  variables  of  the  simulation  such  as  the  total  energy 
and  length  of  each  sequence.  Figure  7  shows  the  string  evolution  under  mesh  refinements. 
The  left  plot  of  Figure  8  shows  the  change  of  total  length  and  total  kinetic  energy  of  the 
string  as  functions  of  time.  The  right  plot  of  the  figure  shows  the  error  of  each  variable  in  the 


sequences.  We  also  integrated  the  norm  of  error  in  both  2D  and  3D  over  time  and  presented 
in  Table  (1)  and  Table  (2).  The  numerical  tests  show  that  the  sequences  are  indeed  hrst 
order  convergent. 

Stress  analysis  on  the  parachute  canopy  and  string  chord  during  the  inflation  process  is 
very  important  to  the  field  test  of  the  parachute.  The  natural  stress  on  each  side  of  the 
triangle  in  the  spring  model  is  the  restoring  force  due  to  the  stretching  of  each  triangle  side. 
Let  Ti,  T2,  ts  be  the  natural  stress  on  side  1,  2, 3  of  a  triangle,  we  have 


Ti  =  k{U-Ll),  z  =  l,2,3, 


(1) 


where  k  is  the  spring  constant,  Lg  is  the  equilibrium  length  of  side  i  and  L*  is  the  stretched 
length  of  the  side  i.  This  natural  stress  can  be  converted  []  into  the  stress  in  Cartesian 
coordinates  on  the  plane  of  the  triangle.  The  Cartesian  stress  is  a  2  x  2  tensor  in  the  plane 
of  the  triangle 


The  conversion  from  natural  stress  to  Cartesian  stress  is  through  a  mapping  matrix,  that  is 


SiCi 

S2C2 

S3C3 


(2) 


where  q  =  cos  Oi  =  dxi/Li  and  s,  =  sin6*i  =  dyi/Li  are  the  cosine  and  sine  functions  of  the 
angle  each  side  with  the  x— axis.  The  stresses  in  two  principal  directions  are  obtained  via 
diagonalization  of  the  stress  tensor,  that  is 

f  CTl  0  ^ 

V  0  a2  ;  ’ 

where  ui  and  (72  are  the  solutions  of  the  characteristic  equation 


(^1,2 


(^xy 


^xy  ^yy  *^1,2 


=  0. 


We  use  the  von  Mises  stress 


'^vm  =  Y  +  cri  -  cric^2  (3) 

to  measure  the  tension  on  the  surface  of  the  fabric.  The  safety  factor  of  the  material  is 
defined  as  the  ratio  of  the  significant  strength  of  the  material  to  the  von  Mises  stress.  When 
this  factor  is  decreased  to  a  critical  value  (which  corresponding  to  high  stress  of  the  fabric 
surface),  it  sends  a  warning  signal  for  the  design  of  the  parachute  canopy.  Figure  9  shows 
the  computed  stress  of  rectangular  and  triangular  membranes  as  we  implemented  during  the 
summer  visit. 

We  had  two  meetings  with  the  parachute  testing  team  led  by  Alec  Dyatt.  The  parachute 
team  made  some  vital  comments  on  the  direction  of  the  parachute  study.  They  suggested  that 
our  parachute  validation  should  include  shape  comparison,  attitude  dependence,  the  varying 
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angle  of  chute  skirt  with  the  airflow  on  the  opening  of  the  parachute,  and  the  quantification  of 
opening  force  on  parachute  payload.  Under  their  suggestion  we  obtained  the  data  of  inflation 
drag  on  the  C-9  parachute  from  Prof.  Jean  Potvin  at  St  Louis  University.  Figure  10  shows 
the  comparison  of  our  front  tracking  simulation  in  comparison  with  Potvin’s  data  on  C-9 
parachute. 


C-9  Canopy  -  Comparison  with  PCPRG  Drops 

{Drops  performed  in  2&04  by  J.  Potvin  andG.  Peek} 


Figure  10;  Evolution  of  drag  during  the  inflation  phase  of  C-9  personnel  parachute.  The  experimental  data 
is  provided  by  Dr.  Jean  Potvin  at  St.  Louis  University. 


4.  Study  of  Airplane  Icing 

In-flight  icing  is  the  result  of  super-cooled  water  on  the  airframe.  It  can  be  in  the 
form  of  cloud  droplets  or  freezing  rain/drizzle.  Aircraft  icing  can  adversely  affect  the  flight 
characteristics  such  as  increase  drag,  decrease  lift,  and  can  cause  control  problems.  In  worst 
case,  it  can  cause  the  airplane  to  crash. 

The  Air  Force  testing  center  has  the  charter  to  test  aircraft  and  weapons  in  all  weather 
conditions  including  in  icing  conditions  and  in  rain.  They  use  the  AIT  (airborne  icing  tanker) 
to  perform  the  artificial  in-flight  icing  and  rain  testing.  These  tests  are  aimed  to  ensure  that 
aircraft  and  weapons  have  capability  to  operate  in  the  climatic  extreme  conditions. 

The  system  consists  nozzles  for  generating  jets  spray  in  the  condition  similar  to  icing  and 
rain  condition  in  flight.  The  testing  is  costly  but  necessary.  Computer  simulation,  although 
cannot  replace  such  tests,  can  provide  useful  information  and  carry  out  certain  virtual  testing 
at  much  lower  cost. 

At  Edwards  AFB,  we  communicated  with  the  base  scientists.  Dr.  Reasor  and  Mr. 
Bhamidipati  about  the  capability  of  front  tracking  code  in  the  study  of  phase  transition 
problem  such  as  deposition,  dissolution,  freezing  and  melting  [4].  Figure  11  shows  examples 
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Figure  11;  Front  tracking  simulation  of  jet  spray,  crystal  formation,  and  collective  motion  of  droplets. 
The  inter-operable  modules  in  the  front  tracking  library  is  an  ideal  computational  tool  for  airplane  icing 
simulation. 

of  front  tracking  modules  for  simulations  of  jet  spray  and  crystal  formation.  The  interop¬ 
erability  of  the  front  tracking  modules  makes  it  capable  for  the  simulation  of  airplane  icing 
test. 
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